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Using Monte Carlo methods and finite-size scaling, we investigate surface crit- 

icality in the O(n) models on the simple-cubic lattice with n = 1, 2, and 3, 

i.e. the Ising, XY, and Heisenberg models. For the critical couplings we find 

K c (n = 2) = 0.4541655(10) and K c (n = 3) = 0.693002(2). We simulate the 

three models with open surfaces and determine the surface magnetic exponents at 

the ordinary transition to be = 0.7374 (15), 0.781 (2), and 0.813 (2) for n = 1, 2, 

and 3, respectively. Then we vary the surface coupling K\ and locate the so-called 

special transition at K c {n = 1) = 0.50214(8) and n c (n = 2) = 0.6222(3), where 

k = Ki/K — 1. The corresponding surface thermal and magnetic exponents are 

y^i = 0.715 (1) and yj^ = 1.636 (1) for the Ising model, and y[f = 0.608 (4) and 

y^l = 1.675 (1) for the XY model. Finite-size corrections with an exponent close 

to —1/2 occur for both models. Also for the Heisenberg model we find substantial 

evidence for the existence of a special surface transition. 

PACS numbers: 05.50.+q, 64.60.Cn, 64.60.Fr, 75.10.Hk 

I. INTRODUCTION 

In the past decades, surface effects near a phase transition have been investigated exten- 
sively, and many results have been obtained by means of mean-field theory, series expansions, 
renormalization and field-theoretic analyses. For reviews, see e.g. Refs. 1,2, and for more 
recent work see Refs. 3,4. In particular, at a second-order phase transition, where long-range 
correlations emerge, surface effects can be significant. The surfaces display critical phenom- 
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ena which differ from the bulk critical behavior; several surface universality classes can exist 
for one bulk universality class. We shall refer to the various types of transitions using the 
terminology of Ref. 1. 

In this work, we investigate surface critical phenomena in three-dimensional O(n) models, 
namely the Ising (n = 1), the XY (n = 2), and the Heisenberg (n = 3) model. The reduced 
Hamiltonian of these models can be written as the sum of two parts: a bulk term proportional 
to the volume of the system and a surface term proportional to the surface area, i.e., 

H/k B T = -K2_^ si- sj - H • s k - K 1 2_^ s p ■ s q - H l ■ ^ s r , (1) 

(ij) k (pq) r 

where the dynamic variable s is a unit vector of n components. The parameters K and 
K\ are the strenghts of the coupling between nearest-neighbor sites in the bulk and on the 
surface layers, respectively, and H and Hi represent the reduced magnetic fields. The first 
two sums in Eq. (1) account for the bulk and the last two sums involve the spins on the 
open surfaces. For ferromagnetic bulk and surface couplings (K > and K 1 > 0), the 
phase transitions are sketched in Fig. 1 for the case of the Ising and the XY model. In the 
high-temperature region, i.e., for bulk coupling K < K c , the bulk is in the paramagnetic 
state, so that the bulk correlation length remains finite. However, a phase transition can still 
occur on the open surface when the surface coupling K\ is sufficiently enhanced. This phase 
transition is referred to as the "surface transition," and is represented by the solid curve in 
Fig. 1. These phase transitions are generally thought to be in the same universality classes 
as the two-dimensional Ising and the XY model, respectively. At the bulk critical point 
K = K c , the line of surface phase transitions terminates at a point (K c , K\ c ). At this point, 
both the surface and the bulk correlation length diverge. Thus, the point (K c , K\ c ) acts 
as a multicritical point, and the phase transition is referred to as the "special transition." 
For Ki < Ki c , the bulk and the surfaces simultaneously undergo a phase transition at 
K = K c . In this case, the critical correlations on the surfaces arise from the diverging bulk 
correlation lengths, and the transition is named the "ordinary transition." The ordinary 
transition remains within the same universality class for a wide range of surface couplings. 
The correlation functions on and near the surface appear to fit universal profiles 5 . The 
transitions at K = K c for K 1 > K lc are referred to as the "extraordinary transitions." 
For the Ising model, since the surfaces are already in the ferromagnetic state for a smaller 
coupling K < K c , no surface transition occurs when the bulk critical line K = K c is crossed. 
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FIG. 1: Sketch of the surface phase transitions of the three-dimensional Ising and XY models 
with ferromagnetic couplings. The vertical axis is the bulk temperature 1/K, and the parameter 
k = {K\ — K)/K in the horizontal axis represents the enhancement of the surface couplings. The 
"surface," the "ordinary," and the "extraordinary" phase transitions are represented by the thick 
solid, the thin solid, and the dashed line, respectively. The lines meet in a point, shown as the 
black circle, which is referred to as the "special" phase transition. 

Nevertheless, owing to the diverging bulk correlation length, the surfaces still display critical 
correlations at K = K c . For the XY model, however, the surface transitions for K < K c 
are Kosterlitz-Thouless-like 6 , i.e., the surfaces do not display long-range order for K < K c , 
in agreement with results of Landau and co-workers 7 . 

For three-dimensional O(n) models with n > 2, which include the Heisenberg model, 
the line of surface transitions for K < K c does not exist; it may thus seem self-evident 
that the special and the extraordinary transitions are also absent. However, this remains 
to be investigated; for instance, in two-dimensional tricritical Potts models, a line of edge 
transitions is absent, but special and extraordinary transitions do exist 8 . Thus, even without 
a line of surface transitions for K < K c , rich surface critical phenomena can still occur in the 
three-dimensional Heisenberg model. For instance, it was reported 9 that at bulk criticality 
K = K c the surface magnetic exponents depend on the ratio K\/K for Ki/K > 2.0. This 
brings up the question whether 
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Additional surface critical phenomena can occur for the Ising model, if the surface and/or 
the bulk couplings are allowed to be antiferromagnetic. Further, one can allow the spins on 
the surface to vanish, such that the surface part of the Hamiltonian in Eq. (1) is described 
by the so-called Blume-Capel model. Such spin-0 states act as annealed vacancies on the 
surfaces. It was observed 10 that, by varying the fugacity of the vacancies, one can reach 
a point where the bulk Ising criticality K = K c joins the line of surface transitions that 
belongs to the universality class of the two-dimensional tricritical Ising model. This point 
was named 10 the "tricritical special" phase transition. In short, for each bulk universality 
class, surface transitions in various surface universality classes can occur, including the 
ordinary, special, and extraordinary transitions at K — K c , and the surface transitions at 
K<K C . 

Apart from the bulk renormalization exponents, additional surface exponents are needed 
to describe the above surface critical behavior. At the ordinary and the extraordinary 
transitions, the surface magnetic scaling field is relevant, while the surface thermal field is 
irrelevant. At the special transition, both the magnetic and the thermal surface fields are 
relevant. 

Since exact information about critical behavior is scarce in three dimensions, deter- 
minations of these surface critical exponents rely on approximations of various kinds. 
These include mean- field theory 1 ' 11 ' 12 ' 13 , series expansions 14 , renormalization group tech- 
nique 2 ' 3,15 ' 16 ' 17 , Monte Carlo simulations 5,7 ' 18 ' 19 ' 20 ' 21 ' 22 , etc. 

The surface critical index f3\ is defined so as to describe the asymptotic scaling behavior 
of the surface magnetization mi as a function of the bulk thermal field t, i.e., mi oc t^ 1 . 
From the scaling relations it follows that this exponent is related to the critical exponents 
as Pi = (d — 1 — Uhi)/yt, where y t and yhi are the bulk thermal and the surface magnetic 
exponent, respectively, and d = 3 is the spatial dimensionality. The mean- field analysis and 
the Gaussian fixed point of the 4 theory yield the magnetic surface index Pi as p[°^ = 1, 
p[ s ^ = 1/2, and pf 1 = 1 respectively for the ordinary, special, and extraordinary transition. 
Many numerical results also exist. For the simple-cubic lattice, the special transition of the 
Ising model was located as k c = 0.5004 (2) 19 ' 20 . Although the values of critical couplings K c 
and Ki c are far from the mean-field predictions, the above result for k c is in agreement with 
the mean- field value k c — 1/2. Further, the surface critical exponents are determined 19 ' 20 ' 21 ' 23 
as = 0.737(5), = 1.62(2), and y^f = 0.94(6). Compared to the Ising model, 
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there are fewer investigations for the three-dimensional XY and the Heisenberg model. 
In particular, to our knowledge, numerical determinations of the special transition and the 
corresponding surface critical exponents have not yet been reported for the XY model. Most 
of the existing results for the Ising, the XY and the Heisenberg model will be tabulated 
below, together with new results of the present work. 

The present work aims to provide an extensive and systematic Monte Carlo investigation 
of the phase transitions of the surfaces of the three-dimensional Ising, XY, and Heisenberg 
models. Compared to numerical investigations one or two decades ago, one has the following 
advantages. Firstly, the bulk critical points of these systems have now been determined accu- 
rately. On the simple cubic lattice, the bulk critical point of the Ising model was determined 
as K c (n = 1) = 0.221 654 55 (3) 24 , with the uncertainty only in the eighth decimal place. The 
bulk transitions of the XY and the Heisenberg model were also determined 14 ' 25 ' 26 ' 27,28 ' 29 ' 30 
to occur at K c = 0.454 167 (4) and 0.693 002 (12), respectively. In the present paper, we also 
simulate these two systems with periodic boundary conditions, and improve the above esti- 
mations as K c (n = 2) = 0.454 1655 (10) and K c (n = 3) = 0.693 002 (2). Secondly, the rapid 
development of computer technology makes it possible to perform extensive computations 
at a limited cost. The present work was performed on 20 PCs; the total computer time is 
in the order of 20 CPU months at a processor speed of 2.5 GHz. 

The organization of the present paper is as follows. Section II reviews the finite-size- 
scaling properties of the systems defined by Eq. (1), with the emphasis on the sampled 
quantities required for the numerical analysis of the simulation data. Section III describes 
the determination of the critical points of the XY and Heisenberg models. Sections IV, 
V, and VI present the Monte Carlo simulations and the results for the Ising, XY, and 
Heisenberg models, respectively. Section VII concludes the paper with a brief discussion. 

II. FINITE-SIZE SCALING AND SAMPLED QUANTITIES 

The total free energy of a system with free surfaces can, in analogy with the Hamiltonian 
in Eq. (1), be expressed as the sum of a bulk and a surface term 1 ' 31 ' 32 : 



F = f h V + hS , 



(2) 
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where / b and f\ are the densities of the bulk and the surface parts of the free energy, 
respectively, and V and S represent the total volume and the surface area, respectively. 
Near criticality, the finite-size scaling behavior of /b and /i is given by the equations 

/ b (t, h, L) = L~ d f hs (tL yt , hV») + / b a(t, h) , (3) 

and 

h{t,hMM,L)=L- {d - 1) h s {tL y \hL y \t 1 L y *\h 1 ^^ . (4) 

The functions /b s and / ba are the singular and the analytical parts of / b ; fi s and /i a similarly 
apply to the surface free-energy density f\. The bulk thermal and magnetic scaling fields 
are represented by t and h, and the surface scaling fields by ti and h\. The associated 
exponents are labeled with corresponding subscripts. As implied by Eq. (3), the leading 
scaling behavior of the bulk does not depend on the presence of free surfaces, although 
physical quantities near the surfaces can be significantly affected. 

On the basis of Eqs. (3) and (4), the scaling behavior of various quantities can be obtained 
as derivatives of / b and / s with respect to the appropriate scaling fields. Details can be found 
in Ref. 1. 

The determination of the bulk critical points used simulations of L x L x L with peri- 
odic boundary conditions in which case fi vanishes. The sampling procedure involved the 
determination of the bulk magnetization density 

N 

rh = N- l Y,s k , (5) 
fc=i 

where N = L 3 . This yielded the averages of the magnetization moments (m ■ m) and 
((m • m) 2 ). The quantity 

Q (K,L) = p%, (6) 
{(m ■ my) 

which is related to the Binder cumulant 33 , converges to a universal value Q at the critical 
point, and was used to determine the critical coupling K c . The finite-size scaling behavior of 
Q can be found by writing the moments of rh in terms of derivatives of the free energy with 
respect to the magnetic field. After application of a scaling transformation, the singular 
powers in Q associated with field derivatives cancel, as do the powers of the nonuniversal 
metric factor relating the physical field and the magnetic scaling field. In the vicinity of 



7 



the critical point one obtains, in terms of the temperature scaling field t and an irrelevant 
temperature-like field u, 

Q(t,u,L) = Q(tL y \uL yi ) + b 2 L^ 2 ^ + b z U^ + ■■■ (7) 

where y x is the leading irrelevant exponent. The correction term with amplitude b 2 is due 
to the analytic contribution to the second moment of rh, and that with amplitude b 3 to 
the second-order dependence of the temperature field on the physical magnetic field. Apart 
from corrections, the temperature field is proportional to K — K c . Eq. (7) will be used in 
Sec. Ill to determine the bulk critical points. 

In order to investigate surface critical behavior, we simulated L x L x L simple-cubic 
lattices with periodic boundary conditions in the xy plane and free boundaries in the z 
direction. First, we sampled the components of the surface magnetization and obtained two 
generalized surface susceptibilities: 

L 2 

Xn = — (m x • m x + m 2 • m 2 ) , and X12 = L 2 (rhi ■ rh 2 ) (8) 

where rh\ and m 2 are the magnetization densities at the free surfaces with z — 1 and z = L, 
respectively. By differentiating the surface free energy with respect to magnetic fields that 
act on either one of the free surfaces, one finds that the singular parts of these surface 
susceptibilities scale as L 2yhl ~ 2 . 

In addition, we computed two surface- surface correlations 

L 

2L' 2 



l/ri ^T2 ^2 ((^.f. 1 ' $x+r,y+r,l + Sx,y,L • S x+r ^ +r ^)) (r — L/2) , (9) 



x,y=l 

and 



x,y=l 

Further, we sampled two ratios of surface magnetization moments: 

(mi -mi) 2 (mi-m 2 ) 2 , , 

Qn = 7T^ and Q 12 = — — — . (11) 

((m 1 -m 1 J z ) {(mi ■ m 2 ) ) 

These quantities are the surface analogs of the bulk ratio Q, cf. Eq. (7), and will be used to 
locate the surface transitions. 
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TABLE I: Description of the simulations of the XY and Heisenberg models. The table lists the 
simulation length in millions of cycles (#MC), and the number of Wolff clusters (#Wc/C) per 
cycle, for each system size L. The data were taken range AK of the coupling K. The val- 
ues shown are those for the XY model; those for the Heisenberg model are approximately the same. 



L 


#MC 


#Wc/C 


AK 


4 


50 


2 





6 


50 


3 


0.016 


8 


50 


4 





10 


20 


5 





12 


20 


6 





14 


20 


7 





16 


80 


8 


0.006 


20 


20 


10 





24 


20 


12 





28 


20 


14 





32 


80 


16 


0.002 


40 


20 


20 





48 


20 


24 





64 


20 


32 





96 


15 


48 





160 


6.7 


80 






III. CRITICAL POINTS OF THE 0(2) AND THE 0(3) MODELS 

The critical point of the Ising model on the simple cubic lattice is already known 24 with 
sufficient accuracy for the present purposes. We therefore restrict ourselves to the XY and 
Heisenberg models. We used the Wolff cluster algorithm 34 ' 35 to simulate these models on 
simple-cubic lattices with periodic boundary conditions. Each cluster is constructed on the 
basis of one component of the spin vectors. The cluster formation process is thus very similar 
to that for the Ising model. Each simulation consists of a large number of cycles, each of 
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which contains several Wolff steps and a data sampling procedure. The Wolff cluster flips 
do not change the absolute values of the spin components. Thus, to satisfy ergodicity, each 
cycle also includes a random rotation of the whole system of spin vectors. We simulated a 
number of L 3 systems whose finite sizes L are listed in Table I, together with the number 
of Wolff clusters per cycle and the total number of cycles per system size. 

Most simulations of the XY model took place at K — 0.454 15, and of the Heisenberg 
model at K — 0.693. Both values are already very close to the final estimates that we 
shall report for the respective critical points. To avoid bias effects associated with short 
binary shift registers 36 ' 37 we took two such shift registers, with lengths equal to the Mersenne 
exponents 127 and 9689, and added the resulting two maximum-length bit sequences modulo 
2. This procedure leads to a sequence whose leading deviation from randomness is a 9- 
bit correlation, which is a considerable improvement in comparison with the usual 3-bit 
correlations 38 . 

The simulations yielded data for the Binder cumulant as described in the preceding 
Section. Expanding Q in Eq. (7) and expressing the temperature deviation from the critical 
point in K — K c , leads to 

Q(K,L) = Q + ai {K-K c )L yt +a 2 {K-K c ) 2 L 2yt + - ■ .+&iL M + b 2 L 3 ~ 2yh +b 3 L y *~ 2yh + • • • (12) 

where Q is a universal constant and the correction term with amplitude b\ is due to the 
irrelevant field. This expression was used to analyze the numerical data for Q(K, L) by 
means of least-squares fits. The exponents were set to the estimates obtained by Guida 
and Zinn- Justin 39 , namely y t = 1.492, y x = —0.789 and yh = 2.482 for the XY model, and 
y t = 1.414, y\ = —0.782 and yh = 2.482 for the Heisenberg model. In order to determine the 
amplitudes ai and a 2 we included some data for relatively small (L = 8, 16 and 32) systems, 
taken at values of K differing up to the order of one percent from K c . Satisfactory fits, as 
judged from the residual x 2 compared with the number of degrees of freedom, were obtained 
including all system sizes down to L = 4 for the XY and L = 6 for the Heisenberg model. 
We found that mixed terms proportional to (K — K c )L Vi+yt were insignificant. 

The results for the critical points are K c = 0.454 1655 (10) for the XY model and K c = 
0.693 002(2) for the Heisenberg model. The universal values of the amplitude ratios are 
Q = 0.8049 (2) for the XY model and Q = 0.8776 (2) for the Heisenberg model. We 
obtained similar results with other types of fits, which involved fewer correction terms and 
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TABLE II: Summary of recent results for the critical coupling K c of the three-dimensional XY 
and Heisenberg models on the simple-cubic lattice with nearest-neighbor interactions. The error 
margin in the last decimal place is shown in parentheses. 



Reference 


model 


Year 


K c 


Janke 40 


0(2) 


1993 


0.45408 (8) 


Adler et al. 28 


0(2) 


1993 


0.45414 (7) 


Gottlob and Hasenbusch 41 


0(2) 


1993 


0.454 20 (2) 


Butera and Comi 42 


0(2) 


1997 


0.45419 (3) 


Ballesteros et. al 30 


0(2) 


1996 


0.454165 (4) 


Cucchieri et. al 29 


0(2) 


2002 


0.454167 (4) 


Present work 


0(2) 


2004 


0.4541655 (10) 


Chen et al. 43 


0(3) 


1993 


0.693 035 (37) 


Holm and Janke 44 


0(3) 


1993 


0.693 (1) 


Butera and Comi 42 


0(3) 


1997 


0.693 05 (4) 


Ballesteros et al. 30 


0(3) 


1996 


0.693 002 (12) 


Present work 


0(3) 


2004 


0.693 002 (2) 



excluded some of the smallest system sizes so as to obtain satisfactory residuals. This 
internal consistency confirms that our error estimates are realistic. The present results and 
some recent values taken from the literature are summarized in Table II. 

IV. ISING MODEL 

Although the three-dimensional Ising model has not been exactly solved, considerable 
information about its critical behavior is available from extensive investigations using various 
kinds of approximations. For a review see, e.g., Ref. 45 . For instance, evidence has been 
found that the Ising model is conformally invariant in three dimensions 23 ' 46 . There is some 
consensus that the values of the bulk thermal and magnetic exponents, y t and yh, are 
1.587 and 2.482, respectively, with uncertainty only in the last decimal place. The bulk 
critical points of a variety of three-dimensional systems with Ising universality have also been 
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obtained 24 ; the bulk transition of the Ising model with nearest-neighbor interactions on the 
simple-cubic lattice was determined as K c = 0.221 654 55 (3). The present work conveniently 
chooses this model so that no further work to determine K c is required. As mentioned earlier, 
periodic boundary conditions are imposed in the xy plane and free boundaries along the z 
direction. 



A. Ordinary phase transition 

Using the Wolff cluster algorithm 34 ' 35 , we simulated the Ising model at bulk criticality, 
with the surface couplings chosen equal to the bulk couplings, i.e., K x — K — K c . The 
system sizes were taken as 16 even values in the range 4 < L < 48. During the Monte 
Carlo simulations, we sampled the surface susceptibilities Xu an d X12, and the correlation 
functions gn and gi 2 . To estimate yj^, the universal surface magnetic exponent of the 
ordinary surface transition, we modeled the Monte Carlo data for the surface susceptibilities 
Xu and X12 by expressions of the form 

Xi(L) = Xa + L^- 2 (bo + hL yi + b 2 L y ^ + b 3 L y * + b 4 L y *) , (13) 

where Xa and the b\ are non-universal and depend on the characteristics of the surface; 
X\ stands for either one of Xn an d X2i- The various parameters in this expression were 
determined by a least-squares fit. We set Xa = to fit Xn- 

Similarly, we fitted data for the correlation functions gn and g± 2 to expressions of the 
form 

9l (L) = L 2y >S- A (b + hL Vi + b 2 L y ^ + b 3 L y * + b 4 L y *) , (14) 

Again, g\ can be either gn or g± 2 ; the no n- universal amplitudes 6; are fitting parameters 
independent of the corresponding amplitudes in Eq. (13), although we use the same symbols. 

The correction terms with amplitudes b X) b 2 , 63, and 64 in Eqs. (13) and (14) account for 
the leading finite-size corrections. The exponent y x = —0.821 (5) 24 is the leading irrelevant 
thermal scaling field in the three-dimensional Ising universality class. Further, since the 
thermal surface scaling field for the ordinary transition is irrelevant, it may also introduce 
finite-size corrections. From a simple scaling argument it can be derived that the value of 
this irrelevant surface exponent is y^ = — l 47 , independent of the spatial dimensionality. 
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0.004 
0.003 
£ 0.002 
0.001 


0.001 0.002 0.003 

L -2.526 

FIG. 2: Surface correlation function g±2 vs. £~ 2 - 526 for the Ising model with e = 0.8. The error 
margins, in this figure as well as in the following ones, are of the same order as the thickness of 
the lines. 

In principle, finite-size corrections from other sources can occur, so that we also include the 
terms with amplitudes b 3 and 6 4 . We simply took y 3 = —2 and y 4 = —3. 

Separate fits of the Xn an d X12 data, employing Eq. (13), yield consistent estimates: 
yJJ = 0.736 (2) and 0.738 (2), respectively. 

Fits of <7n and g 12 yield = 0.737 (2) and 0.736 (2), respectively. A joint fit of both sets 
of susceptibility data, as well as one of both sets of correlation function data, employing a 
single parameter yj^ and independently variable amplitudes, yielded consistent results but 
no significant improvement of the accuracy. 

We also simulated Ising systems in which the surface enhancement is defined as in Ref. 5. 
These systems differ from Eq. (1) as to the couplings between the surface layer and the 
second layer. We thus introduce an enhancement parameter e and define couplings K 1 = e 2 K 
between nearest-neighbor sites on the surface, and couplings K[ = eK between surface sites 
and their nearest neighbors in next layer. Whenever we parametrize the surface enhancement 
by e we refer to the Hamiltonian defined in Ref. 5, which differs from Eq. (1). 

By varying the parameter e, one can move closer to the fixed point for the ordinary phase 
transition so as to reduce the amplitudes of finite-size corrections. Systems with e = 1 
reduce to those described above. In accordance with Ref. 5, in the present work we also 
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TABLE III: Summary of the results for the surface critical exponents in the three-dimensional Ising 
model, as obtained by different techniques. MF: mean-field theory, MC: Monte Carlo simulations, 
FT: field-theoretical methods, CI: Conformal invariance. The MF values of yn and yhi have already 
made use of the mean-field predictions for the bulk thermal and magnetic exponents, which are 
yt = 3/2 and yh = 9/4, respectively. 





ordinary 


special 




Vhi 


Pi 


Vhi ya 


01 


4> 


MF 1 ' 11 


1/2 


1 


5/4 3/4 


1/2 


1/2 


MC 18 


0.72(3) 


0.78 (2) 


1.71(16) 0.94(5) 


0.18(2) 


0.59 (4) 


MC 20 


0.721 (6) 


0.807 (4) 


1.623(3) 


0.2375(15) 




MC 5 


0.740(15) 










MC 21 


0.73(1) 


0.80(1) 








MC+CI 23 


0.737(5) 


0.798(5) 








MC 19 






1.624(8) 0.73(2) 


0.237(5) 


0.461 (15) 


■prp2,15 


0.737 


0.796 


1.583 0.855 


0.263 


0.539 


FT 16 


0.706 


0.816 








FT 17 






1.611 1.08 


0.245 


0.68 


Present 


0.7374(15) 0.796(1) 


1.636(1) 0.715(1 


) 0.229(1) 


0.451(1) 



chose e = 0.9 and 0.8. The analyses of the data for the surface susceptibilities and the 
correlation functions again employ Eqs. (13) and (14); the results for the surface magnetic 
exponents are in agreement with those obtained for the case e = 1. As an illustration, the 
data for g 12 with e = 0.8 are shown versus L 2j/ ii )_4 in Fig. 2, where y$ = 0.737 is taken 
from the fit. 

Finally, a joint fit to the data for xn an d X12 for the three cases e = 1.0, 0.9, and 0.8 
yields y$ = 0.7374 (15); this result is in good agreement with most of the existing results, 
as shown in Table III. 



14 




0.48 0.49 0.5 0.51 0.52 0.53 



FIG. 3: Surface dimensionless ratio Qn vs. surface-coupling enhancement k for the Ising model. 
The data points +, x, □, Q, A, 0, and * represent system sizes L = 21,25,29,33,41,49, and 63, 
respectively. 

B. Special phase transition 

Since it is known that the special transition is located near k = (Ki/K) — 1 = 0.5, the 
simulations were performed with surface enhancements n in the range from 0.46 to 0.54, in 
steps of 0.01. The system sizes assumed 18 values in the range 5 < L < 95. We sampled 
several quantities, including the surface susceptibilities Xn and X12, and the universal ratios 
Qn and Qi 2 . Part of the data for Qn are shown in Fig. 3, in which the clear intersection 
indicates the location of the special transition. As mentioned earlier, when k deviates 
from the finite-size behavior of Qn is governed by the surface thermal exponent yf x . 
We fitted the data for Q n and Q 12 by 

4 4 

Qi(«, L) = Qg + J- a k {n - KffL^ + b i Lm + 

k=l 1=1 

c(« - K^)L y " +yi + n(n - k®) 2 L v " + r Q L ya + 

n(« - K^)L ya + r 2 (n - K^) 2 L ya + r 3 (K - K^) 3 L y « , (15) 

where the terms with amplitude bi account for various finite-size corrections; and again 
the subindex 1 in Qi and Q± c is shorthand for 11 or 12, whichever the case may be. The 
terms with amplitudes r\ (i — 0, • • • ,3) are due to the analytic background. The derivation 
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of Eq. (15) can be found e.g. in Ref. 24. Naturally, we fixed the exponent y 1 = y x = 
—0.821 (5) 24 , the exponent of the leading irrelevant scaling field in the three-dimensional 
Ising model. In principle, additional corrections due to irrelevant scaling fields can be 
induced by the open surfaces, so that we set y 2 = yn as an unknown exponent. In order 
to reduce the residual x 2 without discarding data for many small system sizes, we included 
further finite-size corrections with integer powers y 3 = —2 and y 4 = —3. The term with 
coefficient n reflects the nonlinear dependence of the scaling field on k, and the one with c 
describes the "mixed" effect of the surface thermal field and the irrelevant field. The terms 
with amplitudes r , n, r 2 , and r 3 arise from the analytical part of the free energy, and the 
exponent y a is equal to 2 — 1y h [- As determined later, the surface magnetic exponent at 
the special transition is about y^\ = 1.636 (1), so that we fixed the exponent y a = —1.272. 
The fits of Q u yields Q llc = 0.626 (1), k£° = 0.50214 (8), and = 0.7154 (14); from the 
fit of Q 12 , we obtain Q 12c = 0.2689(1), re£° = 0.50207(8), and y$ = 0.715(4). Next, we 
simultaneously fitted the data for Qu and Q\ 2 by Eq. (15), and obtain = 0.50208 (5), 
and y^f = 0.715 (1). Our estimate = 0.50208 (5) does not agree well with the existing 
results = 0.5004 (2) 19,20 . Further, as expected, k$ does not assume the mean-field value 
1/2. Attempts to determine the unknown exponent yn and its associated amplitude by 
least-square fitting to the Qu and Q12 data were unsuccessful. These corrections, if present, 
do not exceed the detection threshold. We also fitted the data for the surface susceptibilities 
Xn and X12 by 

4 

X i(k, L) = L 2 *" - 2 [a + ^a k {n- k< 8 ) + b 1 L^ + b 2 L m + + 

k=i 

hL y4 + c(k - K^)L y ti +yi + n ( K _ K (*)) 2 L y ti + r L ya + 

ri(« - Kf)L y * + t 2 (k - K^) 2 L ya + t 3 (k - re^) 3 ^" + 

C2i(« - ^)L y ^ +m + c 22 {k - K ( S ))2 L 2^ ) + Kl ] _ (16) 

Again, the correction exponents were taken as y x = —0.821 (5) 24 , y s = —2, and y 4 = —3, 
and the exponent y 2 = yn was left to be fitted. Other than in Eq. (15), we have included 
in Eq. (16) the combined effect of the surface thermal field and the irrelevant field with the 
unknown exponent yn, as described by the mixed terms with amplitudes c 2 \ and c 22 . These 
terms lead to a reduction of the residual x 2 of the fits, but do not significantly modify the 
result for the surface exponent y^\. The surface thermal exponent was fixed at y^f = 0.715 
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FIG. 4: Surface susceptibility xnL^ 1 ' 272 vs. surface-coupling enhancement k for the Ising model. 
The data points +, x, □, 0> A, 0, and * represent system sizes L = 21,25,29,33,41,49, and 63, 
respectively. 

as found above. The fit of xn yields k£° = 0.50209 (9), = 1.636 (1), and y n = -0.52(2). 
The quoted error margins include the uncertainty due to the error in y^f. In this case we 
found clear evidence for corrections, introduced by the surfaces with an exponent yn. It is 
remarkable that such corrections are significant only in combination with /t-dependent terms. 
The data for the surface susceptibility are shown in Fig. 4 as L)^ 1 - 272 , where the 
exponent, which stands for 2 — 2y h {, is chosen such as to suppress the leading L-dependence 
at the special transition. As expected, the data display intersections approaching the special 
transition as determined above. 

V. XY MODEL 

The bulk critical point of the XY model was determined as K c = 0.454 1655 (10) in Sec. 
II, which is of sufficient accuracy to perform the following simulations only at this estimated 
value of K c . 
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A. Ordinary phase transition 

In analogy with the Ising model, we first let the surface couplings K 1 assume the same 
values of the bulk couplings, i.e., K\ — K — K c . The system size took 14 values in the 
range 4 < L < 48. We sampled the surface susceptibilities xn and X12, and the correlation 
functions gn and g i2 , and analyzed the data as we did for the Ising model at the ordinary 
phase transition. For instance, the data for xu and X12 were also fitted by Eq. (13), in which 
the irrelevant exponent was taken as y^ = — 0.789 39 . The estimates of the surface magnetic 
exponent y^i from various quantities agree; the result is y^i = 0.781 (2). 

As a consistency test, in analogy with the Ising model, we also simulated the surface- 
enhanced XY model as defined in Ref. 5, with e = 0.9 and 0.8. As expected, the results for 
these two cases are in good agreement with the above estimate y$ = 0.781 (2). However, 
since the simulations are less extensive in comparison with those for the case e = 1, they do 
not significantly improve the accuracy of y^ . 

B. Special phase transition 

As discussed above, the XY model is a marginal case in the sense that the line of surface 
phase transitions for K < K c is Kosterlitz-Thouless-like. Still, one would expect that, 
for K = K c , the special and the extraordinary surface transitions occur. Therefore, we 
performed simulations at the estimated bulk critical point as given above, and varied the 
surface enhancement from k = 0.48 to k — 0.68. The system sizes took on 19 values in the 
range 5 < L < 95. The sampled quantities include the surface susceptibilities Xu and X12, 
the correlation functions gn and g± 2 , and the dimensionless ratios Qn and Q± 2 . Part of the 
data for Q i2 are shown in Fig. 5, where the intersection clearly indicates that the special 
transition occurs near k c = 0.622. Further, the increase of the slope of Q as a function 
of finite size L strongly suggests that the surface thermal exponent at k c is larger than 0, 
i.e., that the scaling field associated with k is not marginal at the special transition. The 
data for Qn and Q 12 were fitted by Eq. (15), in which the leading irrelevant exponent was 
fixed at y\ = — 0.789 39 and the exponent y 2 = yn was left free. We obtain Q llc = 0.840 (1), 
Q 12c = 0.379(2), k c = 0.6222(3), and ylf = 0.608(4). The fits of Q u and Q 12 do not 
provide clear evidence for the existence of a term with exponent ya. 
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FIG. 5: Surface dimensionless ratio Qu vs. surface-coupling enhancement k for the XY model. 
The data points +, x, □, 0> A, 0, and * represent system sizes L = 17,21,25,33,41,49, and 63, 
respectively. 

TABLE IV: Summary of the results for the surface critical exponents in the three-dimensional XY 
and Heisenberg models. MC: Monte Carlo simulations, SE: series expansions. 





ordinary 


special 




Vhi 


Vhi yn 


mc pry) 7 


0.74 




SE (XY) U 


0.81 




MC (XY) 5 


0.790(15) 




Present (XY) 


0.781 (2) 


1.675(1) 0.608(4) 


MC (Heisenberg) 5 


0.79(2) 




Present (Heisenberg) 


0.813(2) 





We also fitted the surface susceptibilities xn an d X12 by Eq. (16). We obtain the sur- 
face magnetic exponent as yjff 1 = 1.675(1). Further, we find evidence for new finite-size- 
corrections with exponent yn = —0.44 (4), the major contribution to which comes form the 
mixed terms with amplitudes c 2 i and c 2 2 in Eq. (16). Results for the surface exponents are 
summarized in Table IV. 
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C. Extraordinary phase transition 

TABLE V: Monte Carlo data for the second moment of surface magnetization m\ and the dimen- 
sionless ratio Qn for the three-dimensional XY model with enhancement k = 1. 



L 

2 

Qu 


7 9 11 13 17 21 25 
0.5653(1) 0.5293(1) 0.5037(1) 0.4839(1) 0.4561(1) 0.4364(1) 0.4216(1) 
0.96242(6) 0.96580(6) 0.96878(5) 0.97138(4) 0.97543(3) 0.97835(3) 0.98065(3) 


L 

2 

Qn 


33 41 49 63 71 81 95 
0.4004(1) 0.3859(1) 0.3747(1) 0.3601(1) 0.3540(1) 0.3473(1) 0.3397(1) 
0.98381(3) 0.98601(3) 0.98748(3) 0.98927(3) 0.99004(3) 0.99085(3) 0.99169(3) 



Two-dimensional surfaces of the XY model do not display spontaneous long-ranged sur- 
face order for K < K c , but they are in a ferromagnetic state in the low-temperature region 
K > K c . Thus the onset of long-range order on the surface also occurs at K — K c . This 
differs from the Ising model, where a long-range ordered surface exists for K < K c if k > k c . 
We performed simulations at k — 1 for the critical XY model with the system sizes in the 
range 7 < L < 95. We sampled the second moment of the surface magnetization m\ and 
the ratio Qn] the data for these two quantities are shown in Table V. 

In order to analyze the finite-size data in Table V, one first requires the proper scaling 
formulas. For the extraordinary phase transitions in the XY model, there exists some am- 
biguity, because it is not generally clear whether the surfaces undergo a first or a second 
order transition. Nevertheless, in either case, the surfaces should display some critical sin- 
gularities, arising from the diverging bulk correlation length. Thus, we fitted the m\ data 
by 

ml(L) =ml + L~ 2X ^(b + biL Vl + b 2 L 2 ^) . (17) 

If the transition on the surface is first order at K = K c , the analytical contribution, m 2 , 
assumes a nonzero value. First, we set the exponent y\ = y x = — 0.789 39 . Satisfactory fits 
were obtained for all the m\ data in Table V, with the terms m\ and those with bo and b\ 
only. The results are m a = 0.471 (5), X$ = 0.188 (5), b = 0.65 (1), and &i = 0.35 (5). The 
quality of the fit is shown in Fig. 6. Further, we fitted the data for the ratio Q n by 
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FIG. 6: Surface magnetization in terms of the quantity (mi) 

2 _ blL -i-2 vs _ L -2X^ for the XY 

model at k = 1, where the values xj^ = 0.188(5) and b\ = 0.35(5) were obtained from a least- 
squares fit (see text). 

Qn(L) =Q C + hL- 2X ^ + b 2 L- 2X ^ + b 3 L~ 2X ^ + ^ + b 4 L~ 2X ^ + ^ , (18) 

where the irrelevant exponent is fixed at y\ — y\ — — 0.789 39 . The presence of the exponent 
xj^ is due to the nonzero background contribution m a in the second moment of the magne- 
tization m 2 . We obtain the asymptotic value Q c = 0.9998 (4) m 1. From the results for m a 
and Q c , it seems that the surface transition at K = K c and k — 1 is first order. However, 
it seems also possible that the surface magnetization vanishes only very slowly as the sys- 
tem size L increases, such that the line of extraordinary transitions on the surfaces is still 
Kosterlitz-Thouless-like. Thus, we set m a in Eq. (17) to zero, and fitted the unknown pa- 
rameters including both xj^ and y\ to the m\ data. Indeed, we found that our Monte Carlo 
data for m 2 in Table V can be modeled this way, and we obtain bo = 0.40 (1), b\ = 0.703 (6), 
X$ = 0.0325 (30), and Vl = -0.545 (14). This fit is illustrated by Fig. 7. We also fitted the 
Q data by Eq. (18) with y x fixed at -0.545, and the result for Q c is Q c = 0.9982 (15), which 
is also consistent with 1. In short, our numerical evidence for the surface magnetization of 
the three-dimensional XY model is not sufficient to determine whether the line of transitions 
for K = K c and k > k c is first or second order, but settling this matter convincingly would 
require extensive simulations, well beyond the scope of the present investigation. 
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FIG. 7: Surface magnetization in terms of the quantity (mi) 2 — b\L °- 61 vs. L °- 065 for the XY 
model at k = 1. 

VI. HEISENBERG MODEL 

We simulated the three-dimensional Heisenberg model at K\ = K = K c = 0.693 002 (2), 
as determined in Sec. II. The system sizes were taken in the range 4 < L < 64. The data 
for the surface susceptibilities Xn an d X12, taken at k — 0, were fitted by Eq. (13). Using a 
similar procedure as that for the XY model, we obtain = 0.813 (2) for the ordinary phase 
transition. We also determined the bulk susceptibility Xh and the dimensionless ratios Qu 
and Q12 for a range of larger values of the surface enhancement k. The scaled susceptibility 
XbL 3 ~ 2yh is shown in Fig. 8. The intersections near k ~ 0.8 are very suggestive of a special 
transition. The results for Qu, shown in Figs. 9 and 10, display similar behavior. For k < 0.8, 
Qu converges to a universal constant characteristic of the ordinary transition. For /t>0.8 
the data seem to converge to a /t-dependent value. The overall behavior of the results for 
Qu resembles that of the ratio Q for bulk transitions in the Kosterlitz-Thouless universality 
class, as reported for the triangular Ising antiferromagnet with nearest- and next-nearest- 
neighbor interactions 48 . An alternative interpretation would be a special transition with 
a relevant exponent y^f only slightly larger than 0. A convincing numerical test of the 
Kosterlitz-Thouless nature of the special transition would require simulations beyond the 
scope of the present work. 
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FIG. 8: Critical bulk susceptibility Xb of the Heisenberg model vs. surface enhancement k. The 
data shown along the vertical axis are scaled with a size-dependent factor L 3 ~ 2yh where = 2.482 
is the bulk magnetic exponent. The data points +, x, □, O) A, 0, and * represent system sizes 
L = 16, 20, 24, 32, 40, 48, and 64, respectively. According to the theory, the scaled susceptibility 
XbL 3 ~ 2yh converges with increasing size L to a value that may still depend on k. The intersections 
near k = 0.8 suggest the existence of a "special" phase transition. 

VII. DISCUSSION 

We used Monte Carlo techniques and finite-size scaling in order to obtain new and more 
accurate results for the bulk and surface critical parameters of the three-dimensional Ising, 
XY, and Heisenberg models. At the ordinary phase transitions, we determined the surface 
magnetic exponents as yj$(n = 1) = 0.7374(15), yj$(n = 2) = 0.781(2), and j/jj(n = 
3) = 0.813(2). These values are in a satisfactory agreement with earlier results 5 , namely 
y JJ(n = 1) = 0.740 (15), y$(n = 2) = 0.790 (15), and y$(n = 3) = 0.79 (2), as shown in 
Table IV. Since the bulk thermal exponent y t of the O(n) model decreases with increasing n, 
these results suggest that the surface exponent y^} is a decreasing function of y t . The same 
seems to hold true for the two- and three-dimensional Potts models, as may be concluded 
on the basis of the following evidence. In three dimensions, the surface magnetic exponent 
for the q — > and q — > 1 Potts models are y$ = 2 and 1.0246 (6) 49 , respectively. The 
former model is generally referred to as the uniform spanning tree 50 , while the q — > 1 Potts 
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FIG. 9: Surface dimensionless ratio Qn vs. surface-coupling enhancement re for the 0(3) model. 
The data points +, x, □, Q, A, 0, and * represent system sizes L = 8,12,16,20,24,32, and 
40, respectively. For small surface enhancement re < 0.5, the ratio Qn converges with increasing 
L to a nontrivial value near 0.62, just as expected for the ordinary phase transition. For large 
enhancement re > 1, it seems that the asymptotic value Qn(L — > oo) is different from 1, and 
dependent on re. In the intermediate range 0.6 < re < 0.9, the slope of the Qn data lines increases 
with L. The intersections of these lines seem to converge to a value near re = 0.8. This figure bears 
much analogy with that for the bulk ratio Q of transitions in the Kosterlitz-Thouless universality 
class. 

model reduces to the bond percolation model. For the two-dimensional Potts model, from 
the conformal field theory, the exponent y^i is exactly known as y^ = 2 — 3/(3 — yt) 51 , 
which is a decreasing function of the bulk thermal exponent yt- Further, if one applies the 
above expression to the tricritical branch of the Potts model in two dimensions, one obtains 
that the surface magnetic scaling field is irrelevant at the ordinary phase transition. Starting 
from this observation, it was found 8 that rich surface phase transitions can also occur in 
some two-dimensional systems, although their "surfaces" are just one-dimensional edges. 

In the present work, we also located the special transitions of the Ising and the XY model 
on the simple-cubic lattice, and obtained numerical estimates of the corresponding renor- 
malization exponents. While the surface transition of the three-dimensional XY model is 
Kosterlitz-Thouless-like, and the line of surface transitions connects to the special transition 
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FIG. 10: Surface ratio Qn in the range 0.65 < k < 1.1 for the 0(3) model. The data points +, 
x, □, 0> ^ 0, and * represent system sizes L = 8,16,24,32,40,48, and 64, respectively. The 
apparent convergence of the intersections of the Qn data with increasing system size indicates a 
"special" surface transition near k = 0.80, in agreement with the results in Figs. 8 and 9. 

point, our numerical data did not yield evidence for corrections to scaling due to a marginal 
field at the special transition. 

Finally, we note that the surface-critical behavior of the 0(1), 0(2) and 0(3) models is 
rather dissimilar for large surface enhancements. For the 0(1) model, spontaneous surface 
order exists even below the bulk critical coupling K c ; for the 0(2) model it exists for K > K c 
and possibly for K = K c ; and for the 0(3) model only for K > K c . In line with the bulk 
critical singularity, the 0(n) surface critical behavior is thus seen to become less singular 
with increasing n. This is also evident from our analyses of the special transitions, which 
yield relevant exponents y^f for the 0(1) and 0(2) models but allow a marginal exponent 
for the 0(3) model. Since the lower critical dimensionality of the special transition 1 is 3 
for n > 2, it seems plausible that the range k > k c corresponds with a line of fixed points 
and ^-dependent critical surface exponents, in agreement with an analysis of the surface 
magnetization by Krech 9 . Indeed, the data in Figs. 8 and 9 are suggestive of a Kosterlitz- 
Thouless-like scenario involving a nonuniversal range of Q-values such as found earlier in 
the different context of the Ising triangular antiferromagnet 48 . 
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